//
// Created by phk on 24-6-20.
//
#include "cosine_approx.h"
#include <gmpxx.h>
#include <vector>
#include <cmath>
#include <vector>
#include <iostream>
#include <complex>
#include <tuple>

const char *piStr = "3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211"
                    "7067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954"
                    "9303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072"
                    "6024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414"
                    "6951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891"
                    "2279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523"
                    "8467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958"
                    "5371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978"
                    "0499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865"
                    "8753320838142061717766914730359825349042875546873115956286388235378759375195778185778053217122680"
                    "66130019278766111959092164201989";

mpf_class CosineApprox::NewFloat(double x)
{
    mpf_class y;
    y.set_prec(1000); // Set precision to 1000 bits
    y = x;
    return y;
}

mpf_class CosineApprox::BigintCos(const mpf_class &x)
{
    mpf_class tmp;
    int k = 1000; // Number of iterations
    mpf_class t = NewFloat(0.5);
    mpf_class half = t;

    for (int i = 1; i < k - 1; ++i) {
        t *= half;
    }

    mpf_class s = x * t;
    s *= x;
    s *= t;

    mpf_class four = NewFloat(4.0);

    for (int i = 1; i < k; ++i) {
        tmp = four - s;
        s *= tmp;
    }

    mpf_class cosx = s / NewFloat(2.0);
    cosx = NewFloat(1.0) - cosx;

    return cosx;
}

mpf_class CosineApprox::BigintSin(const mpf_class &x)
{
    mpf_class sinx = NewFloat(1.0);
    mpf_class tmp = BigintCos(x);
    tmp *= tmp;
    sinx -= tmp;
    mpf_sqrt(sinx.get_mpf_t(), sinx.get_mpf_t());
    return sinx;
}

int CosineApprox::maxIndex(const std::vector<double> &array)
{
    int maxind = 0;
    double max = array[0];
    for (size_t i = 1; i < array.size(); ++i) {
        if (array[i] > max) {
            maxind = int(i);
            max = array[i];
        }
    }
    return maxind;
}

genDeg CosineApprox::genDegrees(int degree, int K, double dev)
{
    if (dev == 0) {
        throw std::runtime_error("Division by zero is not allowed.");
    }
    int degbdd = degree + 1;
    int totdeg = 2 * K - 1;
    double err = 1.0 / dev;

    std::vector<int> deg(K, 1);

    std::vector<double> bdd(K);
    double temp = 0.0;
    for (int i = 1; i <= (2 * K - 1); ++i) {
        temp -= log2(static_cast<double>(i));
    }
    temp += (2.0 * static_cast<double>(K) - 1) * log2(2.0 * mPI);
    temp += log2(err);

    for (int i = 0; i < K; ++i) {
        bdd[i] = temp;
        for (int j = 1; j <= K - 1 - i; ++j) {
            bdd[i] += log2(static_cast<double>(j) + err);
        }
        for (int j = 1; j <= K - 1 + i; ++j) {
            bdd[i] += log2(static_cast<double>(j) + err);
        }
    }

    int maxiter = 200;
    int iter;

    for (iter = 0; iter < maxiter; ++iter) {
        if (totdeg >= degbdd) {
            break;
        }
        int maxi = maxIndex(bdd);

        if (maxi != 0) {
            if (totdeg + 2 > degbdd) {
                break;
            }

            for (int i = 0; i < K; ++i) {
                bdd[i] -= log2(static_cast<double>(totdeg + 1));
                bdd[i] -= log2(static_cast<double>(totdeg + 2));
                bdd[i] += 2.0 * log2(2.0 * mPI);

                if (i != maxi) {
                    bdd[i] += log2(std::abs(static_cast<double>(i - maxi)) + err);
                    bdd[i] += log2(static_cast<double>(i + maxi) + err);
                } else {
                    bdd[i] += log2(err) - 1.0;
                    bdd[i] += log2(2.0 * static_cast<double>(i) + err);
                }
            }

            totdeg += 2;
        } else {
            bdd[0] -= log2(static_cast<double>(totdeg + 1));
            bdd[0] += log2(err) - 1.0;
            bdd[0] += log2(2.0 * mPI);
            for (int i = 1; i < K; ++i) {
                bdd[i] -= log2(static_cast<double>(totdeg + 1));
                bdd[i] += log2(2.0 * mPI);
                bdd[i] += log2(static_cast<double>(i) + err);
            }

            totdeg++;
        }

        deg[maxi]++;
    }
    return { deg, totdeg };
}

genNode CosineApprox::genNodes(const std::vector<int> &deg, double dev, int totdeg, int K, int scnum)
{
    if (dev == 0) {
        throw std::runtime_error("Division by zero is not allowed.");
    }
    mpf_class PI;
    PI.set_prec(1000);
    PI.set_str(piStr, 10);

    mpf_class scfac = NewFloat(static_cast<double>(1 << scnum));
    if (scfac == 0) {
        throw std::runtime_error("Division by zero is not allowed.");
    }
    mpf_class intersize = NewFloat(1.0 / dev);

    std::vector<mpf_class> z(totdeg);
    int cnt = 0;

    if (deg[0] % 2 != 0) {
        z[cnt] = NewFloat(0);
        cnt++;
    }

    mpf_class tmp;

    for (int i = K - 1; i > 0; --i) {
        for (int j = 1; j <= deg[i]; ++j) {
            tmp = NewFloat(static_cast<double>(2 * j - 1));
            tmp *= PI;
            tmp /= NewFloat(static_cast<double>(2 * deg[i]));
            tmp = BigintCos(tmp);
            tmp *= intersize;

            z[cnt] = NewFloat(static_cast<double>(i));
            z[cnt] += tmp;
            cnt++;

            z[cnt] = NewFloat(static_cast<double>(-i));
            z[cnt] -= tmp;
            cnt++;
        }
    }

    for (int j = 1; j <= deg[0] / 2; ++j) {
        tmp = NewFloat(static_cast<double>(2 * j - 1));
        tmp *= PI;
        tmp /= NewFloat(static_cast<double>(2 * deg[0]));
        tmp = BigintCos(tmp);
        tmp *= intersize;

        z[cnt] = NewFloat(0) + tmp;
        cnt++;

        z[cnt] = NewFloat(0) - tmp;
        cnt++;
    }

    std::vector<mpf_class> d(totdeg);
    for (int i = 0; i < totdeg; ++i) {
        d[i] = NewFloat(2.0);
        d[i] *= PI;

        z[i] -= NewFloat(0.25);
        z[i] /= scfac;

        d[i] *= z[i];
        d[i] = BigintCos(d[i]);
    }

    for (int j = 1; j < totdeg; ++j) {
        for (int l = 0; l < totdeg - j; ++l) {
            d[l] = d[l + 1] - d[l];
            tmp = z[l + j] - z[l];
            if (tmp == 0) {
                throw std::runtime_error("Division by zero is not allowed.");
            }
            d[l] /= tmp;
        }
    }

    totdeg++;

    std::vector<mpf_class> x(totdeg);
    for (int i = 0; i < totdeg; ++i) {
        x[i] = NewFloat(static_cast<double>(K));
        x[i] /= scfac;
        tmp = NewFloat(static_cast<double>(i)) * PI;
        tmp /= NewFloat(static_cast<double>(totdeg - 1));
        x[i] *= BigintCos(tmp);
    }

    std::vector<mpf_class> c(totdeg);
    std::vector<mpf_class> p(totdeg);
    for (int i = 0; i < totdeg; ++i) {
        p[i] = d[0];
        for (int j = 1; j < totdeg - 1; ++j) {
            tmp = x[i] - z[j];
            p[i] *= tmp;
            p[i] += d[j];
        }
    }

    return { x, p, c, totdeg };
}

std::vector<std::complex<double>> CosineApprox::ApproximateCos(int K, int degree, double dev, int scnum)
{
    mpf_class scfac = NewFloat(static_cast<double>(1 << scnum));
    if (scfac == 0) {
        throw std::runtime_error("Division by zero is not allowed.");
    }
    genDeg gendeg;
    genNode gennode;

    std::vector<int> deg;
    std::vector<mpf_class> x;
    std::vector<mpf_class> p;
    std::vector<mpf_class> c;

    int totdeg;

    gendeg = genDegrees(degree, K, dev);

    deg = gendeg.a;
    totdeg = gendeg.b;
    gennode = genNodes(deg, dev, totdeg, K, scnum);

    x = gennode.x;
    p = gennode.p;
    c = gennode.c;
    totdeg = gennode.totdeg;

    mpf_class tmp;
    std::vector<std::vector<mpf_class>> T(totdeg, std::vector<mpf_class>(totdeg));

    for (int i = 0; i < totdeg; ++i) {
        T[i][0] = NewFloat(1.0);
        T[i][1] = x[i];
        tmp = NewFloat(static_cast<double>(K)) / scfac;
        T[i][1] /= tmp;

        for (int j = 2; j < totdeg; ++j) {
            T[i][j] = NewFloat(2.0) * (x[i] / tmp) * T[i][j - 1] - T[i][j - 2];
        }
    }

    mpf_class maxabs;
    int maxindex;

    for (int i = 0; i < totdeg - 1; ++i) {
        maxabs = abs(T[i][i]);
        maxindex = i;

        for (int j = i + 1; j < totdeg; ++j) {
            if (abs(T[j][i]) > maxabs) {
                maxabs = abs(T[j][i]);
                maxindex = j;
            }
        }

        if (i != maxindex) {
            for (int j = i; j < totdeg; ++j) {
                std::swap(T[i][j], T[maxindex][j]);
            }
            std::swap(p[i], p[maxindex]);
        }

        for (int j = i + 1; j < totdeg; ++j) {
            T[i][j] /= T[i][i];
        }

        p[i] /= T[i][i];
        T[i][i] = NewFloat(1.0);

        for (int j = i + 1; j < totdeg; ++j) {
            tmp = T[j][i] * p[i];
            p[j] -= tmp;

            for (int l = i + 1; l < totdeg; ++l) {
                T[j][l] -= T[j][i] * T[i][l];
            }
            T[j][i] = NewFloat(0.0);
        }
    }

    c[totdeg - 1] = p[totdeg - 1];
    for (int i = totdeg - 2; i >= 0; --i) {
        c[i] = p[i];
        for (int j = i + 1; j < totdeg; ++j) {
            c[i] -= T[i][j] * c[j];
        }
    }

    totdeg--;
    std::vector<std::complex<double>> res(totdeg);
    for (int i = 0; i < totdeg; ++i) {
        double real_val;
        real_val = c[i].get_d();
        res[i] = std::complex<double>(real_val, 0);
    }

    return res;
}
